Dynamic nonlinearity in large scale dynamos with shear 
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ABSTRACT 

We supplement the mean field dynamo growth equation with the total magnetic helicity 
evolution equation. This provides an explicitly time dependent model for alpha quenching in 
dynamo theory. For dynamos without shear, this approach accounts for the observed large scale 
field growth and saturation in numerical simulations. After a significant kinematic phase, the 
dynamo is resistively quenched, i.e. the saturation time depends on the microscopic resistivity. 
This is independent of whether or not the turbulent diffusivity is resistively quenched. We 
find that the approach is also successful for dynamos that include shear and exhibit migratory 
waves (cycles). In this case however, whether or not the cycle period remains of the order of the 
dynamical time scale at large magnetic Reynolds numbers does depend how on how the turbulent 
magnetic diffusivity quenches. Since this is unconstrained by magnetic helicity conservation, 
the diffusivity is presently an input parameter. Comparison to current numerical experiments 
suggests a turbulent diffusivity that depends only weakly on the magnetic Reynolds number, but 
higher resolution simulations are needed. 



Subject headings: MHD - turbulence 

1. Introduction 

The large scale magnetic field of the sun and 
other stars is frequently modeled using a£l dy- 
namo theory (Moffatt 1978; Parker 1979; Krause, 
& Radler 1980; Zeldovich, Ruzmaikin, & Sokoloff 
1983). This theory has been successful in repro- 
ducing the cyclic behavior of solar and stellar ac- 
tivity as well as the latitudinal migration of belts 
of magnetic activity. The cyclic behavior results 
mainly from the shear (the f2 effect); see the refer- 
ences above. Shear also helps producing a strong 
toroidal field, while the a effect remains respon- 
sible for regenerating poloidal from toroidal field. 
Dynamos without (or weak) shear can also gener- 
ate large scale fields, but now the a effect also 
regenerates toroidal from poloidal field (a 2 dy- 
namo). Such dynamos are usually nonoscillatory. 
On the other hand, not all aO dynamos are oscil- 
latory: in oblate geometries (accretion discs and 
galaxies), dynamos tend to be nonoscillatory (e.g., 
Covas et al. 1999). In simple cartesian geome- 



try with periodic boundaries a effect dynamos al- 
ways exhibit migratory waves once shear is strong 
enough. 

The viability of an a effect dynamo has been 
controversial primarily because the nonlinear 
backreaction of the growing magnetic field on the 
dynamo coefficients has not been well understood. 
At the center of the debate is how to incorporate 
the backreaction into the a effect. It is often taken 
to be of the form a = qk<j(B), where ckk is the 
kinematic value and 



q = 1 + aB /B, 



cq 



(1) 



is a lorentzian quenching function. Here, 



B 2 

is the mean field, B cq is the equipartition field 
strength and a is a dimensionless parameter. In 
recent years there has been mounting concern that 
the value of a might actually be of the order of 
the magnetic Reynolds number, R m (Vainshtcin 
& Cattaneo 1992; Gruzinov & Diamond 1994; 
1995; 1996; Bhattacharjee & Yuan 1995; Catta- 
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neo & Hughes 1996), rather than of order unity 
(e.g., Riidiger, & Kitchatinov 1993). In stars, 
R m ^ 10^"'^, so quenching would set in for rather 
weak fields, suggesting that dynamo-generated 
fields should be much below the equipartition field 
strength. This is therefore referred to as 'catas- 
trophic' quenching. 

However, quenchings of the form (1) have tradi- 
tionally been obtained under the assumption that 
the system is in a steady state. This leads to incor- 
rect predictions about the evolution of the mean 
field. Since the magnetic helicity evolution equa- 
tion must be taken into account, and because it is 
time dependent, the functional form for a quench- 
ing also becomes time dependent. Loosely speak- 
ing, a helical field is one with a field-aligned cur- 
rent. While this characterizes primarily the cur- 
rent helicity density, the magnetic helicity is really 
a volume integral which can be associated with the 
topological linkage of flux lines. Both magnetic 
helicity and topological linkage are conserved in 
the non-resistive (large magnetic Reynolds num- 
ber) limit. This imposes a crucial constraint on 
the field evolution in general, and the evolution of 
the a effect in particular. 

There have been a number of important at- 
tempts to incorporate dynamical a quenching 
based on magnetic helicity conservation (Kleeorin 
k Ruzmaikin 1982; Zeldovich et al. 1983; Klee- 
orin, Rogachevskii, & Ruzmaikin 1995; see also 
Seehafer 1996; Ji 1999). Recently Field & Black- 
man (2002, hereafter FB02) derived from Pou- 
quet, Frisch and Leorat (1976), a simple two-scale 
approach whose solutions and physical interpreta- 
tion were shown to agree well with simulations of 
Brandenburg (2001a, hereafter B01). 

Let us define some terms: we refer to the proce- 
dure of obtaining a quenching from a time depen- 
dent differential equation derived from magnetic 
helicity conservation as "dynamical" quenching. 
In this case, a is obtained by solving an explicitly 
time dependent equation. On the other hand, we 
refer to the quenching as fixed form or "algebraic" 
if a is expressed fixed function of B, where 
by fixed we mean that the functional form is kept 
fixed in time, and the only time dependence enters 
through B. Although certain algebraic quenching 
expressions can emerge as a useful approximation 
in specific temporal regimes, we shall show below 
that only the dynamical quenching is consistent 



with the magnetic helicity equation, and the func- 
tional form of a(B) must change with time. The 
lorentzian quenching formula (1), when applied for 
all time, would be an example of a fixed form al- 
gebraic quenching prescription. 

B01 performed three-dimensional simulations 
of a zero shear helical dynamo in a periodic 
domain and found that the field attains super- 
cquipartition field strengths. The initial growth 
phase is kinematic, and it is only at later times 
that a saturation phase emerges. The final super- 
equipartition field strength is reached after a resis- 
tive time scale. This behavior can be reproduced 
empirically by an a 2 dynamo where both a and 
turbulent magnetic diffusivity, r/ t , are catastroph- 
ically quenched according to Eq. (1) with a ~ i? m 
(B01). Below we show however that this form of 
the turbulent electromotive force is not univer- 
sal. The dynamical expression is degenerate in 
special cases where the mean field is current-free 
or force-free; it therefore reproduces catastrophi- 
cally quenched behavior in those cases, but not in 
others. 

For an a 2 dynamo in a periodic box, growth of 
the large scale magnetic field energy is associated 
directly with the growth of large scale magnetic 
helicity. Since the total magnetic helicity can only 
change resistively, growth of the large scale mag- 
netic helicity implies significant growth of small 
scale magnetic helicity of the opposite sign. By 
taking a as proportional to the difference between 
kinetic and current helicities [the "relative helic- 
ity" as derived by Pouquet et al. (1976)] and using 
a two-scale approach, FB02 developed a model for 
the nonlinear dynamical quenching of a, whose 
equation for a is formally equivalent to that of 
Kleeorin & Ruzmaikin (1982) and Zeldovich et al. 
(1983). Using a two-scale approach, FB02 showed 
that the growth of the small scale magnetic he- 
licity augments the current helicity contribution 
to a which ultimately quenches it. The coupling 
between the small and large scale magnetic he- 
licity equations in this two-scale time dependent 
dynamical quenching theory predicts that at late 
times, the quenching of a is consistent with Eq. (1) 
with a ~ i? m , but that at early times the dynamo 
proceeds kincmatically, independent of R m . The 
results agree with the simulations of B01 both in 
terms of the time evolution of the large scale field 
energy and the saturation value. 
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An explicitly time dependent quenching for- 
mula has sometimes been used in mean-field mod- 
els, but the main motivation in many instances 
was to study chaotic behavior in stellar dynamos 
(Ruzmaikin 1981; Schmalz & Stix 1991; Feudel, 
Jansen, & Kurths f993; Covas et al. 1997; 1998). 
Some kind of dynamical a quenching, but with 
an explicit time lag, has previously been invoked 
by Yoshimura (1978) in order to reproduce long 
term behavior with multiple periods. This ap- 
proach was however rather ad hoc and not based 
on magnetic helicity conservation. It was only re- 
cently that Kleeorin et al. (1995; 2000) pointed out 
that the catastrophic quenching of Vainshtein & 
Cattaneo (1992) can emerge from the dynamical 
quenching approach under certain circumstances 
(e.g. when the mean field is current-free or force- 
free) . 

The plan of the paper is as follows. In §§2 and 
3 we present the governing equations in a form 
most suitable for our analytic and numerical treat- 
ment. We then consider limiting cases such as 
early and late time evolution, the effective a dur- 
ing the different stages, and the effects of shear 
(§4) . We then consider the full time evolution nu- 
merically with and without shear, and compare 
with fully three-dimensional simulations (§5). Un- 
like dynamos without shear, dynamos with shear 
can exhibit dynamo waves and thus cyclic behav- 
ior. In addition to being a distinction of astrophys- 
ical relevance, we will see that this distinction is 
important in assessing the role of turbulent diffu- 
sion. Finally, we present possible extensions of the 
model (§6) and present our conclusions (§7). 

2. The dynamical equation for a 

We use the magnetic helicity equation for the 
fluctuating field as an auxiliary equation that 
needs to be solved simultaneously with the mean- 
field dynamo equation. For the a 2 dynamo case, 
the dynamical quenching theory derived below is 
similar to that of FB02, but here we generalize the 
approach to the afl dynamo. 

In a closed or periodic domain the magnetic 
helicity, ( A • B) , evolves according to 

A(A-B) =-2 Wo (J-B), (2) 

where A (with B = V x A) is the magnetic vector 
potential, J = VxB/^ is the current density, r\ 



is the microscopic magnetic diffusivity, and angu- 
lar brackets denote volume averages. We split the 
magnetic field into mean and fluctuating compo- 
nents, i.e. B = B + b (and similarly for all other 
quantities). Mean fields are here defined by av- 
eraging over one or two coordinate directions, de- 
pending on whether the mean field is two- or one- 
dimensional; see below. The evolution of the mean 
magnetic vector potential is given by 

qX — _ _ _ 

— = £ + U x B wmJ - V<f>, (3) 



where £ = u x b is the electromotive force result- 
ing from small scale velocity and magnetic fields, 
and <f) is the electrostatic potential of the mean 
field which can be chosen arbitrarily without af- 
fecting the magnetic field and magnetic helicity. 
From Eq. (3) one obtains an evolution equation 
for the magnetic helicity of the mean field, 

^(A-B)=2(£.B)-2 Wo (J-B), (4) 

and an evolution equation for the magnetic helicity 
of the fluctuating field, 

A(a-b) =-2(£.B)-2 Wo (j-b), (5) 

such that the sum of the two equations becomes 
Eq. (2). We note that U does not enter Eqs (4) 
and (5). A remarkable property of Eq. (5) is 
that it contains no triple moments, in contrast 
to the energy equation for the fluctuating mag- 
netic field, for example. This property allows a 
closure whereby Eq. (5) is solved along with the 
mean-field equations to ensure that the magnetic 
helicity equation (2) is satisfied exactly. 

We now discuss the functional form of 5. In 
mean-field electrodynamics one can show that for 
isotropic homogeneous turbulence (Moffatt 1978) 

£ = aB i] t fi Q J, (6) 

where a will be specified below and rjt is the tur- 
bulent diffusivity. The anisotropics induced by the 
generated large scale field are ignored. In Eq. (6) 
we have ignored a possible contribution from the 
cross helicity effect (Yoshizawa & Yokoi 1993); in 
what follows, we assume that the small scale cross 
helicity is always small. We have measured its 
contribution to £ in a simulation with shear and 
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found it to be <~ 1/20 of the contribution from the 
a effect. 

To proceed, we take the nonlinear a of the 
form originally proposed in Pouquet et al. (1976), 
where the residual (sum of kinetic and magnetic) 
isotropic and homogeneous a effect, 

a = a K + a M - (7) 

is given by 

a K = -\t{ui ■ u), a M = +±r(j -b)/p . (8) 

Angular brackets denote volume averages. (This 
implies that these and other turbulent transport 
coefficients are constant in space.) Non-isotropic 
tensorial generalizations to (8) are known (Klce- 
orin & Rogachevskii 1999; Rogachevskii & Klee- 
orin 2001). In Eq. (8), r is the correlation time, 
(lj ■ u) is the small scale kinetic helicity (with 
u> = V x u), (j-b) is the small scale current helicity 
(with j = V x b//xo, where po is the vacuum per- 
meability), and po is the density which is assumed 
constant. 

To account for the magnetic influence on the 
turbulent magnetic diffusivity, we assume the form 

77t =77to. 9 (B), (9) 

where 

mo = H u2 > ( 10 ) 

is the kinematic value of the turbulent magnetic 
diffusivity, and g(H) is a quenching function (spec- 
ified later) normalized such that g(Q) = 1. We use 
Eq. (10) to eliminate r in Eq. (8), so 

a M = 7?toMo(j • b)/B 2 q . (11) 

with B 2 q = /i p (u 2 ). 1 

In isotropic turbulence, the spectra of magnetic 
and current helicities are related to each other by 
a k 2 factor where k is the wavenumber. To a good 
approximation this also applies in real space to the 
helicities at the two scales of the fluctuating and 
mean fields. In particular, we have 

(a • b) - /x (j • b)/fc 2 = aM 5 c 2 q /(77 t0 fc f 2 ), (12) 

1 In principle, the effective correlation times in the expres- 
sions for «]yi and r? t o (tm an d TKi sa y) could be different. 
This would correspond to replacing _Bg q — > (tm/tk)Bcc\ m 
the final expressions involving B cq . 



where we have used Eq. (11) to relate (a • b) to 
aM- Here, kf is the characteristic wavenumber of 
the fluctuating field. After multiplying Eq. (5) by 
rjtakf j i? 2 q we obtain an evolution equation for aM 
(Klecorin & Ruzmaikin 1982; see also Zeldovich et 
al. 1983; and Klecorin et al. 1995) 

daM |2 /(^ B ) , °m\ , 1qA 
~dT = - 2Vwk [-BT- + R-) ' < 13 > 

where we have defined the magnetic Reynolds 
number as 

Rm = Vto/V- (14) 

This result agrees with that in Kleeorin et al. 
(1995) if their characteristic length scale of the 
turbulent motions at the surface, l s , is identified 
with 2ir/k[ and if their parameter /i is identified 
with 87r 2 ?? 2 /((u 2 )Z 2 ). Note that solving Eq. (3) 
or (4) together with (13) is equivalent to solving 
Eq. (3) or (4) with (5); FB02 solved (4) with (5) 
for the a 2 dynamo. 

The quenching function for magnetic diffusiv- 
ity, (/(B), is uncertain. Cattanco & Vainshtcin 
(1991) proposed a catastrophic quenching formula 
for g(B) in two-dimensional turbulence. Gruzi- 
nov & Diamond (1994) confirmed this, but found 
no quenching in the three-dimensional case, i.e. 
g = 1 for all field strengths. This is in qualitative 
agreement with numerical simulations of Nord- 
lund, Galsgaard, & Stein (1994). Kitchatinov, 
Riidiger, & Pipin (1994) as well as Rogachevskii & 
Kleeorin (2001) found that g oc | B | 1 for strong 
fields and independent of R m . Instead of using 
their detailed functional (and tensorial) formula- 
tions, we adopt here a simple fit formula 

g = (l+g\(B)\/B cq )- 1 (easel), (15) 

which was also used in B01 who found g w 16 for 
runs with different values of i? m . This formula 
matches the asymptotic form of Eq. (20) of Ro- 
gachevskii & Kleeorin (2001) with g — w 
2.25 if (b 2 ) w ^ioA)(u 2 ) is assumed (but it varies 
only little with the level of small scale magnetic 
energy: g = 2.78 if (b 2 ) = is assumed, for exam- 
ple) . In the following we allow for different values 
of g, including g = 0. We emphasize that our 
prescription for the quenching of rjt is not dynam- 
ical, because the quenching depends on (b 2 ) which 
does not obey such a stringent conservation law as 
(a - b), as does that governing a. Nevertheless, for 
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fully helical fields, (b 2 ) and (a-b) are proportional 
to each other. This, as well as earlier work by B01 
and FB02 motivates use of the expression 

g = a/a K (case II), (16) 

which will also be considered below for compari- 
son. 

3. The complete set of model equations 

To summarize our approach, the problem con- 
sists of simultaneously solving the two equations 

OB 

— =Vx (UxB + aB- (77 + ^)^0 J), (17) 
da _ 2 ( a(B 2 ) -T] t fi Q (3-B) a - a K \ 

d* - ~ Wf [ w q + ) ' 

(18) 

where r/t depends on B via Eq. (9). The r/to coeffi- 
cient in Eq. (18) is however constant. (In practice 
we continue solving for A instead of B.) We em- 
phasize that B is spatially periodic and a and r) t 
are spatially uniform. The generalization to non- 
periodic B and spatially varying forms of a and 
r]t is not straightforward and will be discussed at 
the end of the paper. 

Shear (which models differential rotation) can 
be implemented in the form U = (0, Sk^ 1 cos k\X, 0), 
where k\ is the minimum wavenumber. In this 
case the mean field is two-dimensional, i.e. B = 
B(x, z, t), and comparison with corresponding tur- 
bulence simulations of Brandenburg, Bigazzi, & 
Subramanian (2001, hereafter BBS) Brandenburg, 
Doblcr, & Subramanian (2002, hereafter BDS) is 
possible. A simpler model, that we will also con- 
sider, is one with linear shear, U = (0, Sx,0). In 
that case the mean field is one-dimensional, i.e. 
B = B(z,t). Corresponding turbulence simula- 
tions can be carried out in the shearing sheet ap- 
proximation, which allows pseudo-periodic bound- 
ary conditions in x. In that case, however, there 
are currently only the more complicated simula- 
tions relevant to accretion discs (Brandenburg et 
al. 1995), so comparison with the present work is 
difficult. 

The dynamo efficiency is determined by the 
usual dynamo parameters 

C a = OK/(»m>fci), C ^ = S /(VTokl), (19) 



where r/xo = f] + ?7to- In addition we have to spec- 
ify R m , fef , and a parameter ef which measures the 
degree to which the small scale field is helical (de- 
fined in the next section). As a non-dimensional 
measure of kf we introduce Kf = kf/k\. Further- 
more, in the absence of a satisfactory theory for 77 
quenching, the parameter g also has to be spec- 
ified. The problem is therefore completely de- 
scribed by the 6 dimensionless parameters C a , Cq , 
R m , kf, e f , and g. 

When comparing with simulations we may use 
the rough estimate R m w u Tms /(r]kf). A crude 
estimate for C a can be obtained using Eqs (8) 
and (10) together with (u> • u) w fcf(u 2 ), where 
kf = Cfkf, so ax ~ kfVto- It is then convenient to 
define the nondimensional effective wavenumber of 
the fluctuating field, kf = kf/k\, so 

C a w ftf/t, (20) 

where we have introduced the correction factor 

^r/To/r/to = l + i? m 1 , (21) 

where 77x0 = Vto + V- This ratio is just above unity 
for R m > 1. 

In addition to one- and two-dimensional mod- 
els, we shall also consider a one-mode reduc- 
tion that reduces the one-dimensional vector equa- 
tion for A (or equivalently for B) to 2 ordinary 
(complex) differential equations. This procedure 
is similar, although more general, than that of 
FB02 where only two ordinary differential equa- 
tions Eq. (4) and (13) were solved for the a 2 max- 
imally helical dynamo. In our present case, non- 
maximally helical dynamos with shear can also be 
studied. 

4. Preliminary considerations 
4.1. Final field strength 

For helical (or partially helical) fields, the re- 
sulting steady-state field strength is determined 
by (J • B) = 0; see Eq. (2). In terms of mean and 
fluctuating fields this means 

(J-B) = -0-b); (22) 

see Eq. (41) of B01. 

In order to connect the current helicities with 
magnetic energies, we can now define the effective 
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wavenumbers for mean and fluctuating fields more 
precise and write, 

fc m = fc m e m = mo (J • B)/(B 2 ), (23) 

fc f -fc f e f - Mo (j-b)/(b 2 ). (24) 

Here, k m and kf are characteristic wavenumbers 
of mean and fluctuating fields, and e m and Cf are 
the fractions to which these fields are helical. In 
the final state, k m will be close to the smallest 
wavenumber in the computational domain, k\. In 
the absence of shear, e m is of order unity, but can 
be less if there is shear or if the boundary condi- 
tions do not permit fully helical large scale fields 
(see below). In the presence of shear, e m turns out 
to be inversely proportional to the magnitude of 
the shear. The value of kf , on the other hand, is 
determined by small scale properties of the turbu- 
lence and is assumed known. 

Both fc m and kf are positive. However, e m can 
be negative which is typically the case when ok < 
0. The sign of Cf is defined such that it agrees with 
the sign of e m , i.e. both change sign simultaneously 
and hence k m kf > 0. In more general situations, 
fc m can be different from k\. Both k m and kf are 
defined more generally via 

fc^=Mo(J-B)/(A-B), (25) 

fc f 2 = Mo(j-b>/<a-b>. (26) 

Using Eqs (23) and (24) together with Eq. (22) we 
have 

fc m (B 2 ) = Mo (J • B) = -noQ ■ b) = ki (b 2 ), (27) 

which generalizes Eq. (46) of B01 to the case with 
fractional hclicities; see also Eq. (79) of BDS. 

Although Eq. (27) may not be precisely sat- 
isfied in simulations, there is evidence that it is 
most nearly obeyed when the kinetic and mag- 
netic Reynolds numbers are large (B01). In the 
presence of shear, the large scale magnetic energy 
may be time dependent, in which case Eq. (27) is 
expected to apply only on the time average. 

Next, we want to express the final steady state 
values of (b 2 ) = b\ n and (B ) = Bg n in terms of 
B 2 q . Using Eqs (11) and (24) we have first of all 

a M = -r/toMb 2 )/^ 2 (28) 



On the other hand, in the steady state we have 
from Eqs (4) and (6) 

a K + a M - VTk m = 0, (29) 

where ryx = f] + ?7t is the total magnetic diffusivity. 
These two relations yield 

bj n = QK - Vrkm Bl n = a K - i] T k m 

In models where rj t is also quenched, both small 
scale and large scale field strengths increase as rj t 
is more quenched. 

We note that both large and small scale field 
saturation strengths are determined by helicity 
considerations. In case I with g — g nn (for the fi- 
nal state), we have (b 2 )/i? 2 q = (C a — gfink m )L/Rf. 
Making use of the estimate (20), we have 

(b 2 )/S c 2 q « 1 - tg^hjnu (31) 

so (b 2 ) k, £? 2 q in the limit of large C a , i.e. large 
Kf, and/or small k m , which is the case for afl dy- 
namos. Regardless of the value of gfi D , we have 
always 

(B 2 )/(b 2 ) - kf/ Km (32) 

in the final state. 

We reiterate that our analysis applies to flows 
with helicity In the nonhelical case, «k — k m — 
kf = 0, so Eq. (30) cannot be used. Nevertheless, 
one would expect a finite value of (b 2 ) because of 
small scale dynamo action. The present approach 
is not really designed to handle this case: if we 
multiply the first equation in (30) by kf, and then 
set it equal to 0, we would just obtain = for 
that equation. 

Even in the fully helical case there can be sub- 
stantial small scale contributions. Closer inspec- 
tion of the runs of B01 reveals, however, that such 
contributions are particularly important only in 
the early kinematic phase of the dynamo. 2 

We should also point out that in the saturated state, fcf 
is expected to show a weak R m dependence: assuming a 
Kolmogorov spectrum for b between kf and the dissipation 
wavenumber fcj with k^/kf ~ R^ 4 one finds kf = k[R^( 4 . 
During the growth phase, however, the spectrum of b rises 
with k either like fc 3 / 2 (Kulsrud & Andersen 1992) or like 
fc x /3 (Brandenburg et al. 1996, B01), so in cither case the 
spectrum of b is peaked near k^ and fcf « fcj may then be 
a better approximation. 
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4.2. Dependence of k m on a and S 

In order to estimate &g n and -Bfi n) it is essential 
to know the value of k m . We show here that for 
aVt dynamos, k m is inversely proportional to the 
ratio of shear to turbulent velocities. The result- 
ing relation turns out to be well confirmed by the 
numerical model solutions below. 

In a numerical model, the value of fc m can 
be easily calculated for a given mean field using 
Eq. (23). However, in order to understand the 
dependence of k m on a and S we first consider 
a one-dimensional model. We consider the case 
Co. S> C a (the aft approximation), so that we 
can neglect the a effect in the equation for the 
toroidal magnetic field B y . (In the opposite case, 
Co. <C C a , we expect k m w hi.) We employ the 
gauge cj) = XJ ■ A (Brandenburg ct al. 1995), so 
Eq. (3) can then be written as (BBS) 

dA/dt = -SA y £ + aA' x y + t/tA", (33) 

where primes denote ^-derivatives, x and y are 
unit vectors in the x and y-directions, respectively, 
and ?7t = rj + n t is the total magnetic diffusivity. 

In a marginally excited one-dimensional model 
with periodic boundaries, the solution consists of 
traveling waves, so all volume averages are inde- 
pendent of time. In particular, (J • B) and (B ) 
are constant during a magnetic cycle. Therefore, 
a and i]t are constant during the saturated state, 
and we can use linear theory to find A in the form 
Ae lfcz_lwt . Real and imaginary parts of the rel- 
evant eigenvalue (corresponding to growing solu- 
tions) are 

w cyc = Reu = rrrl&(C a Ca/2) 1 / 2 , (34) 

A = Imcj = uj cyc — ?7Tfc m j (35) 

where C a = a/r]Tk m and Cq — S/^rk^ are effec- 
tive dynamo numbers based on k m (not fc m ) during 
the saturated state. The corresponding eigenvec- 
tor is 

A=-(l + i)(Q J /2C Q ) 1 / 2 x + y. (36) 

This yields for k m = Rc(J*B)/|B| 2 the result 

fc m - (2C a C n ) 1/2 k m /(C a + C a ). (37) 

Since the saturated state corresponds to the 
marginally excited state, we have A = 0, which 



implies uj cyc — ^T^m and C a Cn = 2, so 

e m = k m /k m = 2/(C a + C n ). (38) 

Given that we have made the afl approximation, 
i.e. C a <C Cn, we have e m = 2/Cq. Since Cq. = 
S/rjTk^, we have 

e m w 2 VT k 2 JS. (39) 

In the marginal state, w C yc = ^T^m, so we can 
write e m = 2u> cyc / S, which is useful for diagnos- 
tic purposes. Another useful diagnostic quantity, 
considered also in BBS, is the ratio of toroidal 
to poloidal field strength, Q = {(B 2 y ) / '{lf x )) 1/2 . 
In the steady state, Q^ 1 = e m /\/2- These re- 
lations between Q^ 1 , e m , and Co will later be 
compared with those obtained from the one and 
two-dimensional models. 

4.3. The force-free degeneracy 

Although a can only be obtained by solving 
an explicitly time dependent differential equation, 
it will be of some interest to estimate the ef- 
fective values of a both during the growth and 
saturated phases, and to assess whether or not 
a is catastrophically quenched. We first clar- 
ify the potentially misleading finding of B01 that 
the simulation results are empirically described 
by a model where both a and rj t are catastrophi- 
cally quenched. We begin by making the assump- 
tion that the time derivative in Eq. (18) can be 
dropped; see Appendix A for details. This leads 

to 

_ a K + flm^o(J-B)AB 2 q 

l + i? m (B 2 )/Se 2 q 

which we refer to as the adiabatic approximation. 
Equation (40) was first derived by Gruzinov & Di- 
amond (1994; 1995), Bhattacharjee & Yuan (1995) 
and Kleeorin et al. (1995). 

The adiabatic approximation can be applied 
to nonoscillatory dynamos near the final steady 
state. For oscillatory dynamos the adiabatic ap- 
proximation is generally invalid, except in the spe- 
cial case where (£-B) is constant in time (the one- 
dimensional afl dynamos considered below are 
such an example). We note that in the adiabatic 
approximation k{ does not enter explicitly. It only 
enters if one were to calculate (b 2 ) for a given 
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solution. We also note that for an imposed uni- 
form magnetic field we have J = 0, in which case 
Eq. (40) predicts a catastrophically quenched a. 
This is in complete agreement with the numeri- 
cal results of Cattaneo & Hughes (1996). One can 
also see how in the fully helical case a appears to 
be quenched only in a non-resistive way. 

In the special case where the large scale field 
is force free, which was the case in B01, we have 
(J • B)B = (B )J and can then write the mean 
turbulent electromotive force, £ = aB — ^to/^oJ, 
with a from Eq. (40) and constant Tfto, as 



a K + i? m 77t^ (J • B)/B| 



l + i? m (B 2 )/Be 2 q 

a K B 



B - ryto^oJ 



r?to Mo J 



l + R m (B 2 )/BL l + i? m (B 2 )/B e 2 



(41) 



The reformulation allowed us to combine the (J-B) 
term in the a expression, together with the con- 
stant rjto term into a quenching expression for 
rjt; cf. the last term in Eq. (41). If one idcnti- 
fies a with Qk/(1 + -Rm(B }/B 2 q ), and rj t with 

»7to/(l+i?m(B )/B% ), then this a is different from 
that consistent with magnetic helicity conserva- 
tion, but in this particular case it yields the same 
electromotive force. This explains the excellent 
agreement between the fully helical simulations of 
B01 and models with catastrophically quenched a 
and r]t. In cases with shear, for example, this de- 
generacy is lifted and then only Eq. (40), or its 
time dependent generalization (18), can be used. 

If one however does assume from the outset that 
i]t oc a, then magnetic helicity conservation and 
Eq. (40) produce a slightly different resistively lim- 
ited quenching for both a and rj t as seen below, 
which is also not a bad fit to B01 (see FB02). 

Before we go on analyzing the early saturation 
phase we discuss in more detail the effective value 
of a in the fully saturated state. 

4.4. The effective a during saturation 

In this section we consider the value of a in the 
final saturated state. We use the subscript 'fin' for 
final and emphasize that the resulting expressions 
are only valid in the steady state, in which case 
Eqs (7), (11), and (22) yield 



Therefore, as pointed out in FB02, the a also ap- 
pears on the right of (40) and we must manipulate 
to solve for a. This is given by 



a 



a K 



1 + -Rmfffin 



l + i? m (.9fin + i?L/^c 2 q )' 



(43) 



where ga n — g(Ba n ) is the fraction by which the 
turbulent magnetic diffusivity is quenched in the 
final state. To determine the quenching of a in 
the steady state, one must specify g(B) and solve 
for a. We consider the two cases of rj t presented 
in §2. 

For case I, we can see immediately from (43) 
that for large R m , a/a K = (1 + g^BlJB^)- 1 , 
so a is not resistively quenched unless g& n itself is 
quenched to resistively small values. For case II, 
g = a/aK, so we must manipulate (43) further and 
find a quadratic equation for a. The appropriate 
solution for R nl » 1 is 



a = 



_ J «k(1- BlJBlJ for B* n < B e 2 q , 



otherwise, 



(44) 



(J ■B) = -(a-a K )B*J Vt0 . 



(42) 



so a oc rjt is quenched resistively all the way to 
zero if (B ) /B% q > 1 , which is usually the case in 
the simulations. 

Using the assumption that a nontrivial station- 
ary state is reached at late times (justified by sim- 
ulations) several important points are revealed by 
the above results. First, the reason a can only 
be weakly quenched for an unquenched % is that 
otherwise the field would eventually decay through 
the action of rjt, precluding a stationary state in 
the first place. Eq. (43) also shows that a is 
quenched more strongly than rj t when rjt is inde- 
pendent of a, whereas if % oc a then both a and 
7] are quenched in tandem. We therefore expect 
a higher saturation value of the field strength for 
case II of §2 as compared to case I with g = 0. 

Finally, in the case rj t oc a (case II) the fact 
that the saturation ratio of the mean field to the 
equipartition value turns out to be 3> a/a^ at 
saturation is important because it suggests that 
the early time growth must have a less resistively 
limited form of a. We discuss this in the next 
section. In this context, again note how one might 
be misled by uniform field simulations designed to 
measure a, such as those of Cattaneo & Hughes 
(1994), in which the mean field cannot grow. 
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4.5. Early time evolution 

During the early growth phase the magnetic he- 
licity varies on time scales shorter than the re- 
sistive time, so the last term in Eq. (18) can be 
neglected and so «m evolves then approximately 
according to 



dOM 

dt 



-2r) W kf(a K + a M - mkn 



(B ) 



(45) 



We use this equation to describe the early kine- 

2 

matic time evolution when (B ), and hence also 
aM, grow exponentially FB02 showed that the 
early time evolution leads to a nearly R m inde- 
pendent growth phase at the end of which a sig- 
nificant large scale field growth occurs. We now 
derive this in our present formalism. 

In the following discussion we restrict ourselves 
to the case where r\ is small. Magnetic helicity 
conservation then requires that 



(A-B) 



(a-b) (for i<i k in), 



(46) 



where the time t = ikin marks the end of the expo- 
nential growth phase (and the 'initial' saturation 
time i sa t used in B01). This time is determined by 
the condition that the term in brackets in Eq. (45) 
becomes significantly reduced, i.e. — ckm becomes 
comparable to ckk — Vtok m - Using (28) and (46) 
together with (23) and (25), we obtain the mean 
squared field strengths of the small and large scale 
fields at the end of the kinematic phase, 



B 2 

-"eq 



B kin 

B 2 

-"eq 



(47) 

respectively, where we have included the extra cor- 
rection factor 



1=1 + R 



(48) 



which becomes important for intermediate values 
of i? m . This correction factor results from restor- 
ing the aM/flm term from Eq. (18) in Eq. (45). 
Not surprisingly, at the end of the kinematic phase 
the small scale magnetic energy is almost the same 
as in the final state; cf. Eq. (30). However, the 
large scale magnetic energy is still by a factor 
k^/k 2 smaller than in the final state (although 
e m may be somewhat different in the two stages). 



This result was also obtained by Subramanian 
(2002) using a similar approach. Using Eqs (20) 
and (47), we can write 



B 2 

-"eq 



km/ Cm 

Zk{/ef 



1 - 



Kf 



(49) 



which shows that Skin can be comparable to and 
even in excess of B oq , especially when e m is small 
(strong shear). 

As emphasized in FB02, for an a 2 dynamo, the 
initial evolution to i?kin is significantly more op- 
timistic an estimate than what could have been 
expected based on lorentzian quenching. In the 
case of an a£l dynamo, fc m <C fc m , so i?kin can be 
correspondingly larger. In fact, for 



/C{ < fcm/fcf, 



(50) 



the large scale field begins to exceed the small 
scale field already during the kinematic growth 
phase. Using the estimate e m ~ 2/Cn (§4.2), and 
since Cq — Cn during the kinematic stage, we see 
that large and small scale fields become compara- 
ble when 6{Cq > 2fcf/fc m . Combining this with 
the condition for a marginally excited dynamo, 
C a C n = 2, we have e { C n > 2(i/n m ) 1 / 2 w 2. 
In simulations of rotating convection, ef « 0.03 
(Brandenburg et al. 1996); assuming that this rel- 
atively low value of ef is generally valid, we have 
Co ^ £f/2 s=s 60 for the condition above which 
the large scale field exceeds the small scale field 
already during the kinematic growth phase. This 
condition is likely to be satisfied both for stellar 
and galactic dynamos. 

During the subsequent resistively limited sat- 
uration phase the energy of the large scale field 
grows first linearly, 

(B 2 ) w^l^ti-tkh) (fori>£ k in), (51) 

and saturates later in a resistively limited fashion; 
see Eq. (45) of B01. 

In the limit of large R m for the maximally hc- 

2 

lical case, (46) implies that (B ) remains of the 
order of B^ in for times ikin < t <C l/r/kf. Us- 
ing (7), (8), and (12) we also find that toward the 
end of the kinematic regime, a is quenched non- 
rcsistively as 



a/a K w 1 - (B )/Bl m (for t w t kh 



(52) 
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where we have assumed ckk 3> ?7to&m- The fact 
that this is independent of i? m and of the choice of 
rjt , contrasts with the a formulae for the late-time 
regime considered in the previous section. This 
highlights the need for a fully time dependent dy- 
namical theory to understand the time dependence 
of a quenching. 

In the following section we present and discuss 
results from simple mean-field models using dy- 
namical quenching. 

5. The full time evolution 

We will first consider a onc-dimcnsional a 2 fl 
dynamo model with constant shear, U = (0, Sx, 0). 
Such a model is sometimes used to model stel- 
lar dynamo waves traveling in the latitudinal di- 
rection (e.g. Robinson & Durney 1982), where 
(x,y, z) are identified with spherical coordinates 
(r, <fi, — 0). In terms of the mean vector potential 
A(z,t), the uncurled mean-field induction equa- 
tion reads (BBS) 

dA/dt = £- SAy* - Wo J, (53) 

where /ioJ = —d 2 A/dz 2 and A z = 0. In con- 
trast to Eq. (33) we do not make the af2 approx- 
imation. Later we also consider two-dimensional 
models which can be compared with simulations. 
For quick parameter surveys, however, the onc- 
dimcnsional models in the one-mode truncation 
are quite useful. 

5.1. The one- mode truncation 

We first consider the one-mode truncation (k = 
km = ki), i.e. we assume A = Ae lkmZ , where A(t) 
is complex, and solve the set of two ordinary dif- 
ferential equations for A x and A y , 

dA/dt = £ - SA y i - Wo J, (54) 

where /xoJ = k m A. The two components of the 
magnetic field are B x = —\k m A y and B y — ik m A x . 
The electromotive force is £ = aB — rjtHo J, where 
a is given by Eq. (7) and ojm is obtained by solving 
Eq. (13) using (£ ■ B) = Re(£ • B), where aster- 
isks denote complex conjugation. For diagnostic 
purposes we also monitor fc m = Re(J* • B)/|B| 2 . 

In Fig. 1 we plot the evolution of (B ) and (b 2 ) 
for R m = 10 4 , C a = 2, fc f = 5, e f = 1. Ini- 
tially, both quantities grow exponentially at the 
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^7 1 ^ 
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Fig. 1. — The upper panel shows the evolution of 
the dimensionless (B ) and (b 2 ) (solid and dashed 
lines, respectively, and both in units of £> 2 q ) in a 
doubly-logarithmic plot for an a 2 dynamo with 
rjt = const (case I). The adiabatic approximation 

2 

gives significantly smaller values of (B ) by the 
end of the kinematic growth phase. The dash- 

2 

dotted line gives (B ) for case II (% oc a); (b 2 ) is 
essentially the same in cases I and II. Note that 
time is displayed on a logarithmic scale to accom- 
modate comfortably both the kinematic phase and 
the final saturation. The lower panel gives the 
same three curves in a semi-logarithmic plot. The 
inset shows the slow and small (~ 1%) adjust- 
ment of (b 2 ) during the early saturation phase. 
R m = 10 4 , C a = 2, C n - 0, Kf - 10, e f = 1, 
.9 = 0. 

rate A = a^ki — f]Tok 2 - This phase stops rather 
abruptly at t^ n = A -1 ln(Skin/-Bi n i), and turns 
then into a resistively limited growth phase. Note, 
however, that already by the end of the kinematic 
growth phase the large scale field is a certain frac- 
tion of the equipartition field strength, which is 
independent of the magnetic Reynolds number. 
The equipartition time can be obtained by set- 
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ting (B ) = bl in in Eq. (51) and using Eq. (47), 
so 

2r?A&(t - ikin) = (1 - k m /k { ) 2 , (55) 

where we have used, for simplicity, e m = ef = 1. 
Otherwise the expression on the right hand side of 
Eq. (55) would be more complicated. The main 
point was to show that the equipartition time is 
still resistively limited, which is consistent with 
Fig. 1. 

In the example shown in Fig. 1, C a = 2, so the 
dynamo is only weakly supercritical and the final 
field strengths are 6| n w 0.1 x B 2 q and (B 2 ) w B 2 q , 
in agreement with Eq. (30). 

We recall that in the fully helical case, because 
of the force-free degeneracy, the adiabatic approx- 
imation coincides with the catastrophic quenching 
hypothesis (§4.3), and it reproduces the final sat- 
uration phase rather well (B01). For larger values 
of i? m , however, the adiabatic approximation gives 

2 

significantly lower values of (B ) by the end of 
the kinematic growth phase (see the dotted line in 
Fig. 1). This difference increases with increasing 
values of i? m . 

In the case of an a 2 dynamo, the overall evo- 
lution of small and large scale magnetic energy 
is similar, except that the large scale field is in 
general not fully helical (e m <C 1), because the 
toroidal magnetic field can be amplified regardless 
of magnetic helicity. It only relies on the pres- 
ence of a small poloidal field which must still be 
regenerated by the a effect. The final field am- 
plitude can be much larger than for the a 2 dy- 
namo. In particular, as pointed out in §4.5, the 
mean field can exceed the small scale field already 
during the entire kinematic growth phase if Cq is 
large enough. This is the case in the example de- 
picted in Fig. 2, so there is no crossing of the two 
curves as in Fig. 1. 

For case II (i] t oc a), the saturation field 
strength of the large scale field is significantly 
enhanced. Also, because rjr is now much lower 
in the saturated state, the value of e m (and hence 
A; m ) is now strongly suppressed. This is consistent 
with Eq. (38). Furthermore, rjt is now suppressed 
down to the microscopic value, so the dynamo 
period has increased by a factor R m . 




20 40 60 80 100 120 

^7 1 ^ 




20 40 60 80 100 120 

^7 1 ^ 



Fig. 2. — The upper panel shows the evolution of 
the dimensionless (B ) and (b 2 ) for an a 2 fl dy- 
namo (solid and dashed lines, respectively). The 

2 

dash-dotted line gives (B ) for case II (r/t oc a); 
(b 2 ) is essentially the same in cases I and II. The 
inset shows the evolution of fe m . Note that in 
case II, k m w 0. The lower panel shows B y (solid 
line) and B x (dashed line). The cycle frequency 
is of order ifrk 2 , but the cycle amplitude adjusts 
on a resistive time scale. R m = 10 2 , C a = 0.1, 
Cq = 200, Kf = 5, ef = 1, g = 0. In case II (dash- 
dotted line) the dynamo period is very long and 
the amplitude much higher. 

5.2. Comparison with B01 

In the simulations of B01, Runs 1-3 had a mag- 
netic Prandtl number of unity (v/rj = 1) and 
R m w u rms /r/kf varied from 2.4 to 18. The run 
with the highest magnetic Reynolds number was 
Run 5 with R m 100, but here v/rj — 100, so the 
hydrodynamic Reynolds number was unity and 
there was basically no turbulent mixing. Shear 
was absent in those runs, so Cq = 0. The other dy- 
namo parameters can be estimated as C a w Kf/b w 
5; see Eq. (20). Assuming g = 0, i.e. 7/t = t]to, 
we have from Eq. (30) for the final mean squared 



field strength 



Run 5 of BO 1 



BlJB 2 C(l = K f /k m -i = 3.6... 4.0, (56) 

for Runs 1-3. Here we have put /t m = 1 for 
the effective nondimcnsional wavenumber of the 
large scale field. The actual values of B^ n /B 2 q 
are somewhat smaller (for Run 3, for example, 
Bfi n /B% q = 3.6 instead of 3.9). The agreement 
with the theoretically expected value is quite rea- 
sonable even with g = 0. Full agreement could 
in principle be achieved with negative values of g 
(g = —0.13 for Run 3, for example). It is more 
likely, however, that the actual value of C a is 
somewhat less than our estimate kf/u. 

2 

The saturation behavior of (B )(t), as seen in 
the simulations of B01, was already well repro- 
duced by the adiabatic approximation (see Fig. 21 
of B01). This is because the values of i? m are still 
too small to be able to see significant differences 
between dynamical quenching and the adiabatic 
approximation. For Run 5 there is however a no- 
ticeable slow-down in the saturation behavior of 
the field at At « 15 if the large scale field is iden- 
tified with the spectral energy at k — k\ . This be- 
havior is well reproduced by the one-dimensional 
model if R m = 50 is chosen. The nominal value 
of R m is actually around 100, but this is proba- 
bly unrealistic and would also give too high values 
of the kinematic growth rate compared with the 
simulation. The slow-down at At w 15, which is 
also seen in the one-dimensional model, is not re- 
produced in the one-mode truncation. The reason 
for this difference lies in the fact that at early and 
intermediate times, k m = 2 prevails, and only at 
later times fc m = 1 becomes dominant; see Fig. 3. 

In order to test the dynamical quenching the- 
ory more quantitatively, it would be useful to pro- 
duce new simulations with smaller scale separa- 
tion, e.g. fcf = 2 . . . 3, and an initial seed magnetic 
field where only one of several possible large scale 
eigenfunctions are present. 

5.3. a 2 fl dynamos: a parameter study 

When shear is included (Co 7^ 0), toroidal 
field can be regenerated solely by the S'-term in 
Eq. (54). When Co 3> C a , the dynamo efficiency 
is governed by the product C Q Co- This is the 
regime where dynamical and fixed-form algebraic 
quenching lead to very different behaviors. 




1 r 




10 20 30 40 50 60 



model: R m = 50, C a = 5, g = 




Fig. 3. — Comparison of the space-time (or butter- 
fly) diagram from Run 5 of B01 with that from the 
one-dimensional model. Dark (light) shades indi- 
cate negative (positive) values. R m = 50, T> = 5, 
.9 = 0. 

In order to study the effect of changing various 
input parameters we begin with Table 1 where we 
show the results for different values of C a . Note 
that only for small values of C a do the results 
for e m and Lo cyc /S agree with the prediction of 
§4.2. This is simply because the aft approxima- 
tion made in §4.2 is only valid for C a <C Co- As 
C a increases, the field strength increases approxi- 
mately as predicted by Eq. (30). 

Increasing the value of i? m has no effect on the 
field geometry and time scales (Q" 1 , e m , and A 
are unaffected); see Table 2. The field strength 
changes only when i? m is close to unity. Once i? m 
is above a certain value, the results are essentially 
independent of i? m . For g > 0, the field strength 
generally increases, as expected, and the cycle fre- 
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Table 1: The effect of changing C a and Co in an a 2 Q dynamo with dynamical a quenching using the 
one-mode approximation. For all runs, i? m = 20, g = 0. 



C a 




shki 


b lnl B lq 


B lJ B cq 


Q- 1 




^cyc/5 


A/5 


remark 


3.0 








0.42 


2.1 


1.00 


1.00 





2.0/Cn 


a 2 dynamo 


3.0 


2 


42 


0.48 


3.4 


0.58 


0.71 


0.36 


1.1 


C a /C n = 0{l) 


3.0 


20 


420 


0.59 


16 


0.11 


0.12 


0.079 


0.27 


larger Cn 


1.0 


20 


420 


0.19 


9.5 


0.071 


0.10 


0.050 


0.12 


smaller C a 


0.3 


20 


420 


0.042 


2.1 


0.071 


0.10 


0.050 


0.045 


smaller C a 


0.1 


20 


420 








0.071 


0.100 


0.050 





marginal case 


0.1 


200 


4200 


0.019 


9.5 


0.007 


0.010 


0.005 


0.012 


same C a C'n as in line 3 



Table 2: Effect of changing i? m for two different 
values of C a in the one-mode approximation. The 
increase of small and large scale field strength for 
small values of i? m is explained by the large values 
of i. For all the models with C a — 1.0 we find 
Q- 1 = 0.021, e m = 0.014, us cyc /S = 0.015, A/5 = 
0.045, whilst for all the models with C a = 0.1 wc 
find Q- 1 = 0.005, e m = 0.007, u> cyc /5 = 0.003, 
A/5 = 0.011. For all runs, Cn = 300, 5 = 0, 
and S/r/k 2 = bR m Cn varies between 6 x 10 2 and 
3 x 10 5 . 

R m C a 6 fin /B oq B^/B^ 



1000 


1.0 


0.19 


15 


100 


1.0 


0.19 


15 


10 


1.0 


0.21 


17 


1 


1.0 


0.38 


30 


100 


0.1 


0.019 


13 


10 


0.1 


0.021 


15 


1 


0.1 


0.037 


28 



Table 3: Effect of changing g in the one-mode ap- 
proximation. For all runs , Rm — 20 , C a — 0.1, 
C n = 200, 5/r/fc 2 = 4200, b 2 jB 2 q = 0.02, 
A/5 = 0.012. 

g B 2 JB 2 q Q- 1 e m c cyc /5 

9J 0l7 0M0 0.0050 
0.3 24 0.0031 0.0043 0.0021 
1.0 70 0.0011 0.0015 0.0010 

quency decreases; see Table 3. One can also verify 
that e m decreases as g is increased. Models where 
rjt oc a (case II) tend to produce long cycle peri- 
ods if the dynamo is sufficiently supercritical; sec 
Table 4. 



Table 4: Effect of changing Cn for case II in the 
one-mode approximation. The critical value for 
dynamo action is Cn = 20. For all runs, R m = 20, 
C a = 0.1, b 2 jB 2 q = 0.02, 



Cn 


B lJ B eq 


Q- 1 


e m 


^cyc/5 


22 


31 


0.0032 


0.0045 


0.0023 


30 


43 


0.0023 


0.0033 


0.0017 


50 


75 


0.0014 


0.0019 


0.0010 


100 


161 


0.0007 


0.0010 


0.0005 



5.4. A two-dimensional a 2 Vl dynamo 

In order to compare with the simulations of 
BBS, it is important to consider the appropriate 
geometry and shear profile. As in BBS we use si- 
nusoidal shear, U = (0, 5/cf 1 cos k\X, 0), and the 
mean field is B = B(x, z, t); cf. §3. The results are 
shown in Table 5. The calculations have been car- 
ried out using a sixth order finite difference scheme 
in space and a third order Runge-Kutta scheme in 
time. 

As in the simulations of BBS, we have chosen 
negative values of «k, but this choice only affects 
the direction of propagation of the dynamo waves. 
There are dynamo waves traveling in the posi- 
tive z-direction at x = ±n and in the negative 
z-direction at x = 0, which is consistent with the 
three-dimensional simulations. These waves are 
best seen in a space-time (or butterfly) diagram; 
see Fig. 4. Note also that there is an initial ad- 
justment time during which the overall magnetic 
energy settles onto its final value (consistent with 
resistively limited saturation) and the cycle period 
increases by a small amount. 

Compared with the one-dimensional model, the 
values of Q^ 1 and e m are about 30% larger in 
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Table 5: Results from the two-dimensional a 2 f2 dynamo with dynamical a quenching. In the last two rows, 
the results from the simulations of BBS and BDS are given for comparison. Models AG2 and perhaps also 
Rl show some tentative agreement with BBS; the corresponding numbers are shown in bold face. Models s3 
and SI give some tentative agreement with BDS (where S/r}k\ = 1000); the corresponding numbers are 
shown in italics. 







C a 


Cn 


g 


S/r)k 2 


bL/B 2 CQ 


BL/B 2 CQ 


Q- 1 


e m 
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the two-dimensional model, but uj cyc is about 3 
times smaller. This may reflect the fact that in 
the present geometry the upward and downward 
traveling dynamo waves can propagate less freely, 
because they are now also coupled in the x direc- 
tion. 

For comparison, in the simulation of BBS, the 
inferred input parameters for modeling purposes 
are fc^ = 2 (the field varies in x and z), S/rjk 2 = 
2000, R m w u rms Mf « 80, C a w Kf = 1...2. 
The resulting non-dimensional output quantities 
are B 2 JB 2 q « 30, Q" 1 « 0.02, e m « 0.11, 
uj cyc /S = 0.005 . . . 0.010, and \/(r)kl) = 30. 

Model Rl gives, within a factor of two, about 
the right saturation field strength, and also the 
values of Q^ 1 , e m and oj cyc /S agree reasonably 
well with the simulations, but the kinematic 
growth rate is too high. Also, the value of R m 
is probably larger in the simulation where we es- 



timated R m Rj 80. In order to have the right 
growth rate, C a has to be lowered. In order to 
match then the right saturation field strength we 
have to have g rs 3. One such case is Model AG2, 
where i? m = 100. Now cycle frequency, growth 
rate, as well as Q -1 agree reasonably well with 
the simulation. 

The simulation of BDS is more resistive, 
S/r)k 2 = 1000 and R m rj 30, but the re- 
sulting field strength is only somewhat smaller, 
BlJBl q rj 20, whereas Q" 1 rj 0.02, e m rj 0.11 
and LOcyc/S = 0.013 . . .0.015 are all somewhat en- 
hanced relative to BBS. Models sl...s3 (where 
g = 1) and SI. . . S3 (where g = 3) arc now ap- 
propriate for comparison, because they all have 
S/rjk 2 = 1000. Model SI with g = 3 gives the 
best agreement for _B^ n , but the cycle frequency is 
too small. For Model s3 with g = 1, uj cyc is about 
right, but now £>| n is too small. 
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Fig. 4. — Space-time (or butterfly) diagram of B y 
for Model SI with dynamical quenching. At x = 
—ir (x = 0) the shear has attained and negative 
(positive) maximum. Dark (light) shades indicate 
negative (positive) values. R m = 30, C a = 0.35, 
Co = 33, m = 5, €{ = 1, g = 3. 




dynamical 
adiab approx 



Fig. 5. — Evolution of the large scale magnetic 
energy for Model AG2 (solid line). The dotted 
line gives the comparison with the corresponding 
adiabatic approximation (see appendix A). R m = 
100, C a = 0.5, Cb = 20, K { = 5, et = 1, g = 3. 



In all models the values of 6| n in table 5 are 
smaller than in the simulations. As discussed in 



§4.1, this is readily explained by the fact that our 
model does not take into account small scale dy- 
namo action resulting from the non-helical com- 
ponent of the flow. 

Comparing the two simulations with different 
values of i? m (BBS and BDS), the cycle frequency 
changes by a factor compatible with the ratio of 
the two magnetic Reynolds numbers. This is not 
well reproduced by a quenching expression for rj t 
that is independent of R m (case I). On the other 
hand, if r\ t oc a (case II), cj cyc becomes far smaller 
than what is seen in the simulations. A possi- 
ble remedy would be to have some intermediate 
quenching expression for r] t . We should bear in 
mind, however, that our current model ignores the 
feedback from the large scale motions. Such feed- 
back is indeed present in the simulations, which 
also show much more chaotic behavior (e.g. Fig. 8 
of BBS) than our model; see Fig. 5. A more realis- 
tic model should therefore allow for more degrees 
of freedom. In particular, the quenching should 
be allowed to be nonuniform in space. This and 
other extensions of the model are discussed in the 
next section. 

6. Possible extensions of the model 

The dynamical quenching model allows us now 
to test a number of additional aspects and prop- 
erties that have been (or can be) seen in direct 
simulations. 

6.1. Cross helicity evolution and large 
scale velocity feedback 

Although we were justified in ignoring the small 
scale cross helicity contribution to £, we found 
from the simulation of BBS that the large scale 
cross helicity is non-negligible. This turned out to 
be the result of the forcing function for the large 
scale velocity and the asymmetry of the large scale 
field with respect to x = 0, and thus with respect 
to the large scale velocity. The correlation between 
the forcing function for the large scale flow and the 
large scale magnetic field serves as a driver in the 
cross helicity evolution equation. In principle, we 
should explicitly couple the equation for the cross 
helicity into the model. Our models with imposed 
shear, however, did not produce the required sym- 
metry breaking that would lead to a significant 
contribution to the large scale cross helicity. This 
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is because we treat the large scale velocity as be- 
ing kinematic. This should be explored further in 
future work by accounting for the dynamical feed- 
back from the large scale motions. 

6.2. Antiquenching 

As long as the feedback from the small scale 
motions onto a and % involve the quantity (a • b) , 
i.e. as long as Eqs (4) and (5) remain fully cou- 
pled, Eq. (2) is guaranteed to be satisfied. Thus, 
the magnetic helicity equation is obeyed regard- 
less of how (a • b) is coupled to a or rj t . It 
may be that under certain conditions, a and T] t 
may even increase with increasing field strength, 
which we refer to as 'anti-quenching'. Branden- 
burg, Saar, & Turpin (1998) used such models to 
explain the increase of relative stellar cycle fre- 
quency with increasing field strength. As an il- 
lustrative example, we have considered the case 
a = «k — cum, i.e. with the opposite sign as in 
Eq. (7), and r] t = r? t0 (l + (b 2 ) 2 /B 2 q ), where (b 2 ) 
is linked to (a • b) via Eqs (24) and (26). Note 
that 77 antiquenching goes with the 4th power, so 
that it eventually dominates over a antiquench- 

2 

ing. The resulting evolution of (B ) shows the 
expected resistively limited saturation phase. Al- 
ternatively, if a is made to increase with increasing 
small scale field strength, the resulting large scale 
field can saturate faster than usual. A similar be- 
havior can be modeled by choosing g < 0, which 
also speeds up the initial build-up of large scale 
magnetic energy. This is possible, and consistent 
with the magnetic helicity equation, because the 

2 

initial build-up of (B ) happens in this case simul- 
taneously with a sharp burst in (b 2 ) such that the 
sum of (A • B) and (a • b) is still approximately 
constant. 

6.3. Magnetic buoyancy 

In solar and galactic dynamo theory the pos- 
sibility of rising magnetic flux tubes contributing 
to the a effect has been discussed (Leighton 1969; 
Ferriz-Mas, Schmitt, & Schussler 1994; Hanasz & 
Lesch 1997; Brandenburg & Schmitt 1998; Moss, 
Shukurov, & Sokoloff 1999; Thelen 2000; Spruit 
2002). For the sun, the idea is that flux tubes 
emerge from the toroidal magnetic field belt at the 
bottom of the convection zone and become twisted 
by the Coriolis force. We point out that Eq. (18) 



does already capture part of this effect. If there is 
a strong partly buoyant magnetic field at the bot- 
tom of the convection zone, it would contribute to 
(J • B) and therefore, through Eq. (18), to a. Of 
course, this effect cannot constitute a dynamo on 
its own as there is no source of magnetic energy. 
However, in conjunction with shear, from which 
energy can be tapped, this effect could lead to dy- 
namo action. Modeling this in the framework of 
dynamical quenching would be a suitable way to 
include the effects of magnetic buoyancy such that 
magnetic helicity conservation is obeyed. 

6.4. Oscillatory imposed fields 

If there is a uniformly imposed magnetic field 
that is oscillatory in time, Eq. (18) would predict 
that a is also oscillatory. If the oscillation fre- 
quency is high enough, the adiabatic approxima- 
tion breaks down. One might wonder whether this 
would be a way to test explicitly the dynamical a 
quenching concept numerically. However, it turns 
out that the resulting averaged a is even smaller 
on average than what is predicted based on the 
adiabatic approximation. Thus, although dynam- 
ical quenching enhances dynamo generation in the 
case of self-generated fields, it actually lowers a in 
the presence of imposed oscillatory fields. 

6.5. Selective decay 

In is interesting to note that Eq. (18) can also 
be applied to the case of decaying magnetic fields. 
In that case, it predicts reduced turbulent decay 
if (J • B) 7^ 0. Thus, accordingly, a fully helical 
magnetic field should only decay at the resistive 
rate, whereas a non-helical magnetic field would 
decay at the turbulent diffusive rate. The slow 
decay of helical fields is well known and leads to 
the so-called Taylor states where magnetic helicity 
is maximized and magnetic energy minimized (e.g. 
Montgomery, Turner, & Vahala 1978). 

6.6. Hyperdiffusion 

Numerical experiments allow one to understand 
the physics described by the equations by modi- 
fying certain terms. Particularly enlightening has 
been the use of hyperresistivity (or hyperdiffusion) 
by which the ordinary diffusion operator, rf7 2 , is 
simply replaced by 772 V 4 . The diffusion at small 
scales is usually fixed by the mesh resolution and 
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kept unchanged, but with hyperdiffusion the diffu- 
sion at large scales can be decreased substantially. 
This method is frequently used in turbulence re- 
search, but the effects on helical dynamos arc quite 
striking: the saturation field strength is consider- 
ably enhanced and the saturation phase prolonged 
(Brandenburg & Sarson 2002). Our model repro- 
duces these features if ijkf is replaced by r\ik\, 
i]J is replaced by — 772 V 2 J, and R m = f]t/(r]2kf) 
is used. The simulations of Brandenburg & Sar- 
son (2002) showed (for the helical dynamo with- 
out shear) that the large scale field behavior de- 
pends on the diffusion at the scale of the large scale 
field itself and not, as one might naively expect, 
on the diffusion at small scales. This behavior 
is clearly reproduced by the dynamical quenching 
model: reducing the microscopic r\ in the mean 
field equation (and not in the dynamical quench- 
ing equation) increases saturation time and sat- 
uration value as expected. This can be taken as 
additional validation of the dynamical quenching 
model. 

6.7. Losses of small scale field 

Open boundaries may provide a means of shed- 
ding magnetic hclicity and thereby alleviating the 
magnetic helicity constraint (Blackman & Field 
2000; Klceorin ct al. 2000; 2002). Numerical sim- 
ulations have shown, however, that when no ad- 
ditional boundary physics is imposed to transport 
preferentially quantities of a particular scale, most 
of the magnetic helicity is lost by the large scale 
field (Brandenburg & Dobler 2001). In this case, 
the growth of the large scale field cannot be ac- 
celerated. In order to check whether accelerated 
growth of large scale fields is at least in principle 
possible we have modeled the preferential shed- 
ding of small scale fields in two different ways, 
both with similar results. Adding an overall loss 
term of the form — ay\_jT\ 0%% on the right hand side 
of Eq. (18) leads to substantial increase of the 
large scale field [note that this is distinct from the 
— a M /T in Kleeorin ct al. (2000), which is just the 
same as our second term in Eq. (13)]. Likewise, 
setting (b 2 ) (and hence au) to zero in sporadic in- 
tervals accelerates the growth phase and enhances 
the saturation value. Similar results have mean- 
while also been obtained by restarting Run 3 of 
B01 after sporadically removing magnetic field at 
and below the forcing scale (see Fig. 14 of BDS). 



This confirms an important prediction from the 
dynamical quenching model. 

6.8. Generalization to nonuniform a 

Our approach is based on magnetic helicity 
which is a volume integral. However, in astrophys- 
ical bodies kinetic and magnetic helicities are not 
spatially constant and change sign at the equator. 
Generalizing Eq. (13) to the case of space depen- 
dent qk and q;m seems at first glance straight- 
forward: omit the angular brackets and replace 
d/dt by d/dt 7 as was done already in the early 
work of Kleeorin & Ruzmaikin (1982). It may also 
be necessary to include a local phenomenological 
magnetic helicity flux transport term, for exam- 
ple of the form 77 q V 2 «m (in addition to whatever 
global flux terms may be present, e.g. Blackman & 
Field 2000; Kleeorin ct al. 2002). In the presence 
of large scale (meridional) flows, it may further- 
more be appropriate to use the advective deriva- 
tive, D/Di = d/dt + XJ- V. However, the magnetic 
helicity density, a - b, is not gauge- invariant and it 
is no longer strictly related to j b locally. The hope 
would be that the generalization outlined above 
may still be useful as an approximation. 

We have performed calculations with = 
axosinfciz and C a = «ko/( ? 7t^i) = 5. The re- 
sulting large scale field strength is approximately 
equal to B cq , and depends only weakly on R m . 
This is consistent with Kleeorin et al. (2002). Sim- 
ulations with the same sinusoidal a profile (Bran- 
denburg 2001b) have shown that the resulting 
large scale field varies mostly in the x direction, 
which is incompatible with the present model. 
Also the resulting field strength was actually sig- 
nificantly below B cq . 

In the presence of open boundaries the present 
model predicts large scale field strengths that de- 
crease inversely proportional with R m . This is 
even steeper that what was found in the simu- 
lations (Brandenburg & Dobler 2001). Thus, in 
its present form the dynamical quenching model 
does not reproduce satisfactorily the numerical re- 
sults when a varies in space. However, with the 
help of simulations it should be possible to identify 
which of the steps in the derivation of dynamical a 
quenching are no longer satisfied, and hence what 
the cause of the problem is. 
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7. Conclusions 

The magnetic helicity evolution equation is a 
constraint that must be satisfied by any dynamo 
theory. When we apply this in the mean field for- 
malism with the prescription that the a effect is 
proportional to the difference between kinetic and 
current helicities, dynamical a quenching emerges 
as the only theoretically consistent approach to 
a quenching. This is supported by comparisons 
with numerical simulations of dynamos with and 
without shear. Fixed-form, algebraic quenching 
prescriptions may apply in a specific parameter 
regime (e.g. the saturated phase), but are invalid 
for earlier times and are inconsistent with results 
from time dependent analyses. Only dynamical 
quenching has predictive power. 

A key result from dynamical quenching is that 
near-equipartition large scale field strengths are 
reached independently of the magnetic Reynolds 
number by the end of the kinematic phase. Final 
saturation is only reached on a resistively limited 
(R m dependent) time scale but with a saturation 
value independent of R m and equal to B^ n /B^ = 
kf/k m ; see Eqs (30) and (32). 

Although magnetic helicity conservation pro- 
vides a basis for a dynamical quenching of a, the 
form of r]t must be prescribed at present. We 
have shown that current simulations of a 2 (shear- 
free) dynamos constrain the dynamical quench- 
ing of £ ■ B which is a combination of a and rjt, 
but they do not separately constrain r/t- On the 
other hand, cycle periods, emerging only in dy- 
namos with shear, can. At present, the shear dy- 
namo simulations are best described by a dynam- 
ical quenching theory in which r] t is only weakly 
dependent on the magnetic field; g w 3 in Eq. (15). 
Higher resolution simulations are needed to verify 
this. 

The two-scale dynamical nonlinear quenching 
approach based on magnetic helicity conservation 
discussed herein constitutes an improvement over 
fixed-form algebraic quenching approaches. Nev- 
ertheless, there are aspects of high- Reynolds num- 
ber afl dynamos that may require the theory to be 
augmented. For example, we have only considered 
a spatially uniform a coefficient. Allowing for spa- 
tial gradients in a will introduce local helicity flux 
terms that are important for astrophysical bodies 
where a changes sign across the equator. In ad- 



dition, helicity flux across global boundaries was 
also ignored in our calculation, though we know 
that real systems have boundaries. The associated 
boundary magnetic helicity flow may be important 
in coupling the dynamo growth to magnetic helic- 
ity evolution. 
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A. The adiabatic approximation 



We discuss here the justification for when the time derivative in the dynamical quenching expression can 
be neglected (the adiabatic approximation). 



In the steady state, cum and (B ) are given by Eqs (29) and (30), respectively. Linearizing Eqs (17) 
and (18) about this state yields 



i dq 

2 dt 



-(a K /fcm - Vt)kj 
-(a K /k m - ?7T)fcm 



Vkf 




(Al) 



where q is the state vector for the logarithmic departure from equilibrium. For excited 

solutions, the terms in the first column of the matrix in (Al) are usually positive, even for a£l dynamos. 
Inspecting the diagonal terms shows that near the saturated state, om is adjusting rapidly on a dynamical 

2 

time scale whilst (B ) is marginal and adjusts only indirectly (via «m) on a resistive time scale. We can 
therefore use the adiabatic elimination principle (e.g. Haken 1983) to remove the explicit time dependence 
of «m by replacing Eq. (18) by 







(a(B 2 ) 



?7t^o(J • B)j /B 2 q + a M /R B 



(A2) 



Substituting a M = a- (*K and solving for a leads to Eq. (40). 

The adiabatic approximation corresponds to the limit in which memory effects become negligible. This 
is best seen by considering the integral form of Eq. (18), 



a 



2rjkf 



Jo 



)[a K + i?m%A*o 



(j-b; 

B 2 cq 



with the Green's function 



G(t,t') = exp 




(A3) 



(A4) 



As long as the field is weak, the width of the Green's function is the resistive time scale, but when (B ) / B 2 q 
is of order unity the large i? m factor becomes important and the width of the Green's function reduces to a 
dynamical time scale. In that case, the B dependent terms in brackets can be pulled out of the integrals in 
Eqs (A3) and (A4), in which case Eq. (40) is recovered. 
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